#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 13 19:08:28 2022
@author: zq
"""

import numpy as np 
from scipy import optimize as op 
from scipy.sparse import coo_matrix
import matplotlib.pyplot as plt

def sprase_rand(m, n, p):
    '''
    Imports
    -------
    import numpy as np
    from scipy.sparse import coo_matrix
    
    Parameters
    ----------
    m : int
        Rows of generate sprase-matrix.
    n : int
        Columns of generate sprase-matrix.
    p : float
        Proportion of non-zero elements in sparse-matrix.

    Returns
    -------
    scipy.sparse.coo.coo_matrix
        Returns a random generated sprase-matrix(m x n) M,
        try print(M) to see it or use M.toarray() for trans 
        M to an array. (Gaussian distribution)
        
    Version: 1.0 writen by z.q.feng @2022.03.13
    '''
    if type(m) != int or type(n) != int:
        raise TypeError('Rows(m) and Columns(n) must be an interger!')
    if p <= 0 or p > 1:
        raise ValueError('p must in (0, 1] !')
    # Counts of non-zero elements in sparse-matrix
    count = int(m * n * p)
    # Indexs of non-zero elements in sparse-matrix
    rows = np.random.randint(0, m, count)
    cols = np.random.randint(0, n, count)
    # Values of non-zero elements in sparse-matrix
    # (Gaussian distribution)
    data = np.random.randn(len(rows))
    return coo_matrix((data, (rows, cols)), shape=(m, n))

def OLE_linprog(A, b):
    '''
    Ax = b (Given A & b, try to derive x)
    
    Parameters
    ----------
    A : array like.
    b : array like.

    Returns
    -------
    x : array
        Minimal L1 norm solution of the system of equations.
        
    Reference
    ---------
    YAO Jain-kang. An Algorithm for Minimizing l1-Norm to Overdetermined Linear Eguations[J]. 
    JIANGXI SCE7NICE, 2007, 25(1): 1-4.(Available at: 
    http://d.g.wanfangdata.com.cn/Periodical_jxkx200701002.aspx)
        
    Version: 1.0 writen by z.q.feng @2022.03.13
    '''
    A, b = np.array(A), np.array(b)
    if np.size(A, 0) < np.size(A, 1):
        raise ValueError('Matrix A rows must greater than columns!')
    m, n = A.shape
    # Trans A into two matrix(n x n and (m - n) x n)
    A1, A2 = A[:n, :], A[n:, :]
    # inverse of A1 
    A1_ = np.linalg.inv(A1)
    # c_ij = A2 * inv(A1)
    c = np.dot(A2, A1_)
    # r(n+1:m) = A2*inv(A1)*r(1:n) + d
    d = np.dot(c, b[:n]) - b[n:]
    # Basic-Pursuit, target function = sum(u, v)
    t = np.ones([2 * m, 1])
    # Aeq_ = [c I(m-n)]
    Aeq_ = np.hstack([-c, np.eye(m - n, m - n)])
    # Aeq[u v] = Aeq_ * (u - v)
    Aeq = np.hstack([Aeq_, -Aeq_])
    # u, v > 0
    bounds = [(0, None) for i in range(2 * m)]
    # r0 = [u; v]
    r0 = op.linprog(t, A_eq = Aeq, b_eq = d, bounds = bounds,
                    method='revised simplex')['x']
    # Minimal L1-norm residual vector, r = u - v
    r = np.array([r0[:m] - r0[m:]])
    # Solving compatible linear equations Ax = b + r
    # Generalized inverse solution
    x = np.linalg.pinv(A).dot(b + r.T)
    return x

if __name__ == '__main__':
    m, n = 256, 128
    # 256x128矩阵，每个元素服从Gauss随机分布
    A = np.random.randn(m, n)
    # 128x1矩阵，每个元素服从Gauss随机分布
    # b = A * u, 则 Au - b = 0，u 即为我们要求的精确解
    u = np.random.randn(n, 1)
    b = np.dot(A, u)
    alpha = OLE_linprog(A, b)
    # 恢复残差 
    print('\n恢复残差：', np.linalg.norm(alpha - u, 2))
    # 绘图
    plt.figure(figsize = (8, 6))
    # 绘制原信号
    plt.plot(u, 'b', linewidth = 2, label = 'Original')
    # 绘制恢复信号
    plt.plot(alpha, 'r.-', linewidth = 2, label = 'Recovery')
    plt.grid() 
    plt.legend()
    plt.show()